%% Plot comparison
close all;
clear;
load('/nfs/satmag_work/mnair/projects/compare_hdgm/2012_splines/ABG_1995_2012.ncspline.mat')
abg1 = load('/nfs/satmag_work/mnair/projects/compare_hdgm/2014/2014ObsOutput/ABG2013.txt');
abg2 = load('/nfs/satmag_work/mnair/projects/compare_hdgm/2014/2014ObsOutput/ABG2014.txt');
time_hdgm = datenum(floor(abg1(:,1)),0,0) + (abg1(:,1)-floor(abg1(:,1)))*365;
time_array_min = (datenum(2009,1,1,0,0,30): 1: datenum(2012,12,31,23,59,30))';
f = figure(1);
set(f,'Position',[ 26           1        1875        1094]);
subplot(311);
[breaks,coefs,l,k,d] = unmkpp(sp_x);
pp2 = mkpp(breaks,repmat(k-1:-1:1,d*l,1).*coefs(:,1:k-1),d);
%plot(time_array_min, ppval(time_array_min,sp_x));
ax = plotyy(time_array_min, ppval(time_array_min,sp_x),time_array_min,ppval(pp2,time_array_min));
%set(gca,'FontSize',16);
datetick(ax(1),'x','keeplimits');
set(ax(2),'XTick',[]);
hold;
plot(time_hdgm, abg1(:,2),'r', 'LineWidth' , 2);
plot(time_hdgm, abg2(:,2),'k','LineWidth' , 2);
ylabel('X(nT)');
title('ABG');
datetick('x','keeplimits');
subplot(312);
[breaks,coefs,l,k,d] = unmkpp(sp_y);
pp2 = mkpp(breaks,repmat(k-1:-1:1,d*l,1).*coefs(:,1:k-1),d);
ax = plotyy(time_array_min, ppval(time_array_min,sp_y),time_array_min,ppval(pp2,time_array_min));
datetick(ax(1),'x','keeplimits');
set(ax(2),'XTick',[]);

plot(time_array_min, ppval(time_array_min,sp_y));
set(gca,'FontSize',16);
hold;
plot(time_hdgm, abg1(:,3),'r','LineWidth' , 2);
plot(time_hdgm, abg2(:,3),'k','LineWidth' , 2);
ylabel('Y(nT)');
datetick('x','keeplimits');
subplot(313);
[breaks,coefs,l,k,d] = unmkpp(sp_z);
pp2 = mkpp(breaks,repmat(k-1:-1:1,d*l,1).*coefs(:,1:k-1),d);
ax = plotyy(time_array_min, ppval(time_array_min,sp_z),time_array_min,ppval(pp2,time_array_min));
datetick(ax(1),'x','keeplimits');
set(ax(2),'XTick',[]);


plot(time_array_min, ppval(time_array_min,sp_z));
set(gca,'FontSize',16);
hold;
plot(time_hdgm, abg1(:,4),'r','LineWidth' , 2);
plot(time_hdgm, abg2(:,4),'k','LineWidth' , 2);
legend('observed','HDGM2012','HDGM2013');
datetick('x','keeplimits');
ylabel('Z(nT)');
pause;
% saveas(gcf,'/nfs/satmag_work/mnair/projects/compare_hdgm/2014/figures/ABG' ,'eps2');
% saveas(gcf,'/nfs/satmag_work/mnair/projects/compare_hdgm/2014/figures/ABG','fig');


%% Plot derivatives

close all;
clear;

% selected_obs = ['AAA';	'ABG';	'ABK';	'ASP';	'BDV';	'BEL';	'BFO';	...
%     'BLC';	'BMT';	'BOX';	'CBB';	'CLF';	'CTA';	'CZT';	'DRV';	...
%     'EBR';	'EYR';	'FCC';	'FUR';	'GNA';	'HLP';	'HRN';	'KAK';	...
%     'KDU';	'KNY';	'LRM';	'LVV';	'MAW';	'MBO';	'MCQ';	'MEA';	...
%     'MMB';	'NGK';	'NUR';	'NVS';	'OTT';	'PAF';	'RES';	'SOD';	...
%     'STP';	'STJ';	'THY';	'UPS';	'VIC'];
% 

% selected_obs = ['BMT';'TUC';'ABG';'ASP';'HON';'NGK'];

S = dir('/nfs/satmag_work/mnair/projects/compare_hdgm/2014/2014ObsOutput/*.txt')
for i = 1:length(S),
filestring(i,:) = S(i).name(1:3);
end;

[data_obs, ia, ib ] = unique(filestring,'rows');

S1 = dir('/nfs/satmag_work/mnair/projects/compare_hdgm/2012_splines/*.mat')

for i = 1:length(S1),
filestring1(i,:) = S1(i).name(1:3);
end;
[spline_obs, ia, ib ] = unique(filestring1,'rows');

[selected_obs, ia, ib]= intersect(data_obs, spline_obs, 'rows');



for i = 1:length(selected_obs),
%for i = 1:1,


load(['/nfs/satmag_work/mnair/projects/compare_hdgm/2012_splines/' sprintf('%s',selected_obs(i,:)) '_1995_2012.ncspline.mat']);
abg1 = load(['/nfs/satmag_work/mnair/projects/compare_hdgm/2014/2014ObsOutput/' sprintf('%s',selected_obs(i,:)) '2014.txt']);
abg2 = load(['/nfs/satmag_work/mnair/projects/compare_hdgm/2014/2014ObsOutput/' sprintf('%s',selected_obs(i,:)) '2013.txt']);
abg3 = load(['/nfs/satmag_work/mnair/projects/compare_hdgm/2014/2014ObsOutput_v3/' sprintf('%s',selected_obs(i,:)) '2014.txt']);

% get latitude and longitude of the station
ncfname = ['/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/secular_removed/' sprintf('%s',selected_obs(i,:)) '_1995_2012_Secular_Removed.nc'];
ncid = netcdf.open(ncfname,'NOWRITE');
obj = read_geomag_netcdf_metadata(ncid);
netcdf.close(ncid);

gespatial_lat(i) = obj.geospatial_lat;
gespatial_lon(i) = obj.geospatial_lon;



time_hdgm = datenum(floor(abg1(:,1)),0,0) + (abg1(:,1)-floor(abg1(:,1)))*365;

annual_time_1 = datenum(2001:2013,1,1);
annual_time_2 = datenum(2000:2012,1,1);
actual_sec_time = datenum(2000:2012, 7 , 2);

sec_x_hdgm = interp1(time_hdgm, abg1(:,2), annual_time_1, 'linear' ) - interp1(time_hdgm, abg1(:,2), annual_time_2 );
sec_y_hdgm = interp1(time_hdgm, abg1(:,3), annual_time_1, 'linear' ) - interp1(time_hdgm, abg1(:,3), annual_time_2 );
sec_z_hdgm = interp1(time_hdgm, abg1(:,4), annual_time_1, 'linear' ) - interp1(time_hdgm, abg1(:,4), annual_time_2 );

sec_x_hdgm_v3 = interp1(time_hdgm, abg3(:,2), annual_time_1, 'linear' ) - interp1(time_hdgm, abg3(:,2), annual_time_2 );
sec_y_hdgm_v3 = interp1(time_hdgm, abg3(:,3), annual_time_1, 'linear' ) - interp1(time_hdgm, abg3(:,3), annual_time_2 );
sec_z_hdgm_v3= interp1(time_hdgm, abg3(:,4), annual_time_1, 'linear' ) - interp1(time_hdgm, abg3(:,4), annual_time_2 );


sec_x_hdgm_2013 = interp1(time_hdgm, abg2(:,2), annual_time_1, 'linear' ) - interp1(time_hdgm, abg2(:,2), annual_time_2 );
sec_y_hdgm_2013 = interp1(time_hdgm, abg2(:,3), annual_time_1, 'linear' ) - interp1(time_hdgm, abg2(:,3), annual_time_2 );
sec_z_hdgm_2013 = interp1(time_hdgm, abg2(:,4), annual_time_1, 'linear' ) - interp1(time_hdgm, abg2(:,4), annual_time_2 );


sec_x_obs = interp1(time_hdgm, ppval(sp_x,time_hdgm), annual_time_1, 'linear' ) - interp1(time_hdgm, ppval(sp_x,time_hdgm), annual_time_2 );
sec_y_obs = interp1(time_hdgm, ppval(sp_y,time_hdgm), annual_time_1, 'linear' ) - interp1(time_hdgm, ppval(sp_y,time_hdgm), annual_time_2 );
sec_z_obs = interp1(time_hdgm, ppval(sp_z,time_hdgm), annual_time_1, 'linear' ) - interp1(time_hdgm, ppval(sp_z,time_hdgm), annual_time_2 );

difference_x(i,:) =  sec_x_obs - sec_x_hdgm_v3;
difference_y(i,:) =  sec_y_obs - sec_y_hdgm_v3;
difference_z(i,:) =  sec_z_obs - sec_z_hdgm_v3;

f = figure(1);
set(f,'Position',[ 26           1        1875        1094]);
subplot(311);

plot(actual_sec_time, sec_x_obs,'LineWidth',2);
set(gca,'FontSize',16);
hold on;
plot(actual_sec_time, sec_x_hdgm,'r','LineWidth',2);
plot(actual_sec_time, sec_x_hdgm_v3,'r--','LineWidth',2);
plot(actual_sec_time, sec_x_hdgm_2013,'k','LineWidth',2);
h=legend('OBS','HDGM14 V1','HDGM14 V2','HDGM13');
set(h,'FontSize',8);


datetick('x','keeplimits');
ylabel('X(nT/yr)');
title(sprintf('%s',selected_obs(i,:)));

subplot(312);
plot(actual_sec_time, sec_y_obs,'LineWidth',2);
set(gca,'FontSize',16);
hold on;
plot(actual_sec_time, sec_y_hdgm,'r','LineWidth',2);
plot(actual_sec_time, sec_y_hdgm_v3,'r--','LineWidth',2);

plot(actual_sec_time, sec_y_hdgm_2013,'k','LineWidth',2);

h=legend('OBS','HDGM14 V1','HDGM14 V2','HDGM13');
set(h,'FontSize',8);


datetick('x','keeplimits');
ylabel('Y(nT/yr)');

subplot(313);
plot(actual_sec_time, sec_z_obs,'LineWidth',2);
set(gca,'FontSize',16);
hold on;
plot(actual_sec_time, sec_z_hdgm,'r','LineWidth',2);
plot(actual_sec_time, sec_z_hdgm_v3,'r--','LineWidth',2);

plot(actual_sec_time, sec_z_hdgm_2013,'k','LineWidth',2);

h=legend('OBS','HDGM14 V1','HDGM14 V2','HDGM13');
set(h,'FontSize',8);


datetick('x','keeplimits');
ylabel('Z(nT/yr)');

set(f,'PaperPositionMode','auto');
saveas(f,['/nfs/satmag_work/mnair/projects/compare_hdgm/2014/' sprintf('%s',selected_obs(i,:)) '_XYZ'],'tif');


% saveas(gcf,'/nfs/satmag_work/mnair/projects/compare_hdgm2013/ABG_1stdifference','fig');
% 
% fid = fopen('/nfs/satmag_work/mnair/projects/compare_hdgm2013/ABG_Spline.txt','wt');
% fprintf(fid,'Date     X(nT)     Y(nT)   Z(nT)    dX/dt    dY/dt    dZ/dt\n');
% 
% fclose(fid);

%pause;

clear abg1 abg2 abg3 sp_x sp_y sp_z;
close all;

end;
             
% for i = 1:length(time_hdgm),
%     
%     fprintf(fid,'%7.2f %8.1f %8.1f %8.1f %8.5f %8.5f %8.5f\n',...
%         abg1(i,1),...
%         ppval(time_hdgm(i),sp_x),...
%         ppval(time_hdgm(i),sp_y),...
%         ppval(time_hdgm(i),sp_z),...
%         ppval(ppx,time_hdgm(i)),...
%         ppval(ppy,time_hdgm(i)),...
%         ppval(ppz,time_hdgm(i)));
% end;
%fclose(fid);

%close all;
%

%% Plots maps of diferences

% set LVV and DRV as one since there are issues with fitting

difference_z(15,:) = 1;
difference_y(15,:) = 1;
difference_z(15,:) = 1;

difference_x(27,:) = 1;
difference_y(27,:) = 1;
difference_z(27,:) = 1;

difference_x(28,:) = 1;
difference_y(28,:) = 1;
difference_z(28,:) = 1;

%% 2011.5

subplot(311);

worldmap('world')
showaxes('off');
plotm(lat,long)
for i = 1:43,
plotm(double(gespatial_lat(i)),double(gespatial_lon(i)),'r.','markersize',ceil(abs((difference_x(i,12))))*2);
end;
title('2011.5 Observed X - HDGM (nT/yr)');

subplot(312);

worldmap('world')
showaxes('off');

plotm(lat,long)
for i = 1:43,
plotm(double(gespatial_lat(i)),double(gespatial_lon(i)),'r.','markersize',ceil(abs(difference_y(i,12)))*2);
end;
title('2011.5 Observed Y - HDGM (nT/yr)');

subplot(313);
showaxes('off');

worldmap('world')

plotm(lat,long)
for i = 1:43,
plotm(double(gespatial_lat(i)),double(gespatial_lon(i)),'r.','markersize',ceil(abs(difference_z(i,12)))*2);
end;
title('2011.5 Observed Z - HDGM (nT/yr)');


%% 2010.5

f = figure(1);
set(f,'Position',[  267   105   562   990]);
subplot(311);

worldmap('world')
showaxes('off');
plotm(lat,long)
for i = 1:43,
plotm(double(gespatial_lat(i)),double(gespatial_lon(i)),'r.','markersize',ceil(abs((difference_x(i,11))))*2);
end;
title('2010.5 Observed X - HDGM (nT/yr)');

subplot(312);

worldmap('world')
showaxes('off');

plotm(lat,long)
for i = 1:43,
plotm(double(gespatial_lat(i)),double(gespatial_lon(i)),'r.','markersize',ceil(abs(difference_y(i,11)))*2);
end;
title('2010.5 Observed Y - HDGM (nT/yr)');

subplot(313);
showaxes('off');

worldmap('world')

plotm(lat,long)
for i = 1:43,
plotm(double(gespatial_lat(i)),double(gespatial_lon(i)),'r.','markersize',ceil(abs(difference_z(i,11)))*2);
end;
title('2010.5 Observed Z - HDGM (nT/yr)');


%% Plot comparison

selected_obs = ['AAA';	'ABG';	'ABK';	'ASP';	'BDV';	'BEL';	'BFO';	...
    'BLC';	'BMT';	'BOX';	'CBB';	'CLF';	'CTA';	'CZT';	'DRV';	...
    'EBR';	'EYR';	'FCC';	'FUR';	'GNA';	'HLP';	'HRN';	'KAK';	...
    'KDU';	'KNY';	'LRM';	'LVV';	'MAW';	'MBO';	'MCQ';	'MEA';	...
    'MMB';	'NGK';	'NUR';	'NVS';	'OTT';	'PAF';	'RES';	'SOD';	...
    'SPT';	'STJ';	'THY';	'UPS';	'VIC'];

fid = fopen('/nfs/satmag_work/mnair/projects/compare_hdgm/2014/hdgm_difference.txt','wt');
localdirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/absolute/';


for i = 1:length(selected_obs),
  
    


load(['/nfs/satmag_work/mnair/projects/compare_hdgm/' sprintf('%s',selected_obs(i,:)) '_1995_2012.ncspline.mat']);
abg1 = load(['/nfs/satmag_work/mnair/projects/compare_hdgm/2013/2013ObsOutput/' sprintf('%s',selected_obs(i,:)) '2013.txt']);
abg2 = load(['/nfs/satmag_work/mnair/projects/compare_hdgm/2013/2012OBShdgmOutput/' sprintf('%s',selected_obs(i,:)) '2012.txt']);
load(['/nfs/satmag_work/mnair/projects/compare_hdgm/2013/' sprintf('%s',selected_obs(i,:)) '_1995_2012.ncoriginal_data.mat']);
%time_hdgm = datenum(floor(abg1(:,1)),0,0) + (abg1(:,1) - floor(abg1(:,1)))*365;
time_hdgm = datenum(abg1(:,1),0,0);
time_array_min = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2012,12,31,23,59,30))';

L = isnan(x_data);

time_array = (datenum(2000,1,1,0,0,30): 1: datenum(2012,12,31,23,59,30))';

L =  time_array >= time_array_min(find(~L,1,'first')) & time_array <= time_array_min(find(~L,1,'last'));
time_array(~L) = [];

% f = figure(1);
% set(f,'Position',[ 26           1        1875        1094]);
% 
% subplot(311);
% [breaks,coefs,l,k,d] = unmkpp(sp_x);
% pp2 = mkpp(breaks,repmat(k-1:-1:1,d*l,1).*coefs(:,1:k-1),d);
% plot(time_array_min(2629441:end),x_data(2629441:end),'color',0.8*[1 1 1]);
% hold on
% plot(time_array, ppval(time_array,sp_x));
% set(gca,'FontSize',16);
% datetick(gca,'x','keeplimits');
% %set(ax(2),'XTick',[]);
% hold on;
% plot(time_hdgm, abg1(:,2),'r', 'LineWidth' , 2);
% ylabel('X(nT)');
% title(selected_obs(i,:));
% datetick('x','keeplimits');
% 
% subplot(312);
% [breaks,coefs,l,k,d] = unmkpp(sp_y);
% pp2 = mkpp(breaks,repmat(k-1:-1:1,d*l,1).*coefs(:,1:k-1),d);
% %ax = plotyy(time_array, ppval(time_array,sp_y),time_array,ppval(pp2,time_array));
% plot(time_array_min(2629441:end),y_data(2629441:end),'color',0.8*[1 1 1]);
% hold on;
% plot(time_array, ppval(time_array,sp_y));
% 
% datetick(gca,'x','keeplimits');
% %set(ax(2),'XTick',[]);
% set(gca,'FontSize',16);
% plot(time_hdgm, abg1(:,3),'r','LineWidth' , 2);
% ylabel('Y(nT)');
% datetick('x','keeplimits');
% 
% subplot(313);
% [breaks,coefs,l,k,d] = unmkpp(sp_z);
% pp2 = mkpp(breaks,repmat(k-1:-1:1,d*l,1).*coefs(:,1:k-1),d);
% %ax = plotyy(time_array, ppval(time_array,sp_z),time_array,ppval(pp2,time_array));
% plot(time_array_min(2629441:end),z_data(2629441:end),'color',0.8*[1 1 1]);
% hold on;
% plot(time_array, ppval(time_array,sp_z));
% %set(ax(2),'XTick',[]);
% set(gca,'FontSize',16);
% plot(time_hdgm, abg1(:,4),'r','LineWidth' , 2);
% legend('observed','spline','HDGM2013');
% datetick('x','keeplimits');
% ylabel('Z(nT)');
% 
% 
% saveas(f,['/nfs/satmag_work/mnair/projects/compare_hdgm2013/' selected_obs(i,:) '_XYZ'],'tif');
%saveas(f,'/nfs/satmag_work/mnair/projects/compare_hdgm2013/','fig');

f = figure(2);
set(f,'Position',[ 26           1        1875        1094]);
subplot(211);
f_data = sqrt(x_data.^2 + y_data.^2 + z_data.^2);
f_hdgm2013 = abg1(:,5);
f_hdgm2012 = abg2(:,5);
L = time_hdgm < time_array(1) | time_hdgm > time_array(end);
f_spline = sqrt(ppval(time_hdgm,sp_x).^2 + ppval(time_hdgm,sp_y).^2 + ppval(time_hdgm,sp_z).^2);

plot(time_array_min(2629441:end),f_data(2629441:end),'color',0.8*[1 1 1]);
hold on;
plot(time_hdgm,f_spline,'k');
set(gca,'FontSize',16);
datetick(gca,'x','keeplimits');
%set(ax(2),'XTick',[]);
hold on;
plot(time_hdgm, f_hdgm2013,'b', 'LineWidth' , 2);
plot(time_hdgm, f_hdgm2012,'r', 'LineWidth' , 2);
ylabel('F(nT)');
title(selected_obs(i,:));
datetick('x','keeplimits');
legend('observed','spline','HDGM2013','HDGM2012');



subplot(212);

plot(time_hdgm(~L),f_hdgm2013(~L)-f_spline(~L),'b','LineWidth' , 2);

time_data = time_hdgm(~L);
hdgm2013_resid = f_hdgm2013(~L)-f_spline(~L);
hdgm2012_resid = f_hdgm2012(~L)-f_spline(~L);

L = time_data >= datenum(2001,1,1) & time_data <= datenum(2010,1,1);
median_2001_2010 = median(hdgm2013_resid(L));
nsamples_2001_2010 = sum(L);

median_2001_2010_hdgm2012 = median(hdgm2012_resid(L));
nsamples_2001_2010_hdg2012 = sum(L);


L = time_data >= datenum(2012,1,1) & time_data <= datenum(2013,1,1);
median_2012_2013 = median(hdgm2013_resid(L));
nsamples_2012_2013 = sum(L);

median_2012_2013_hdgm2012 = median(hdgm2012_resid(L));
nsamples_2012_2013_hdgm2012 = sum(L);


L = time_data >= datenum(2011,1,1) & time_data <= datenum(2012,1,1);
median_2011 = median(hdgm2013_resid(L));
nsamples_2011 = sum(L);

median_2011_hdgm2012 = median(hdgm2012_resid(L));
nsamples_2011_hdg2012 = sum(L);

 set(gca,'FontSize',16);
 ylabel('F(nT)');
 datetick('x','keeplimits');
 hold on;
 plot(time_hdgm(~L),f_hdgm2012(~L)-f_spline(~L),'r','LineWidth' , 2);

 legend('Obs-hdgm2013','Obs-hdgm2012');
 title(sprintf('%s %4.0f %4.0f %4.0f %4.0f %4.0f %4.0f',selected_obs(i,:),median_2001_2010,...
     median_2012_2013,median_2012_2013_hdgm2012, median_2011, median_2011_hdgm2012 ));

nc_fname = [localdirectory selected_obs(i,:) '_1995_2012.nc'];
obj = read_netCDF_header(nc_fname);

fprintf(fid,'%s %6.2f %6.2f %6.2f %d %6.2f %d %6.2f %d %6.2f %d %6.2f %d %6.2f %d\n', ...
selected_obs(i,:), ...
obj.geospatial_lat, ...
obj.geospatial_lon, ...
median_2001_2010,...
nsamples_2001_2010, ...
median_2012_2013, ...
nsamples_2012_2013, ...
median_2001_2010_hdgm2012, ...
nsamples_2001_2010_hdg2012, ...
median_2012_2013_hdgm2012, ...
nsamples_2012_2013_hdgm2012, ...
median_2011, ...
nsamples_2011, ...
median_2011_hdgm2012, ...
nsamples_2011_hdg2012 ...
);

set(f,'PaperPositionMode','auto');
saveas(f,['/nfs/satmag_work/mnair/projects/compare_hdgm2013/' sprintf('%s',selected_obs(i,:)) '_F'],'tif');
% saveas(f,['/nfs/satmag_work/mnair/projects/compare_hdgm2013/' sprintf('%s',selected_obs(i,:)) '_F'],'fig');


close all;
end;

fclose all;
